python计算莫兰指数(Moran‘s I)并绘制地区热力图

2024-07-16 18:00| 来源: 网络整理| 查看: 265

文章目录程序简介程序/数据集下载代码分析

程序简介

计算省的pm2.5平均值作为观测矩阵，省会的距离的倒数作为空间权重矩阵，计算全局莫兰指数为0.49，显著性检验p值为3.75>1.96，得出中国地区的pm2.5存在空间正相关输入：中国各地区pm2.5值输出：莫兰散点图、莫兰指数、莫兰检验数、地区热力图 莫兰指数(Moran’s I)是空间自相关系数的一种，其值分布在[－1,1]，用于判别空间是否存在自相关。比如：一个小区内有钱人会集聚住在豪宅别墅区，平民们(比如我- -)则会集聚住在普通的楼房区，这里的个人财富就是一种观测值，呈现高高聚集，低低聚集的现象。

程序/数据集下载

点击进入下载地址

代码分析

导入模块，路径

# -*- coding: utf-8 -*- from Module.CalPosition import calDistance from Module.MoranI import moranI import pandas as pd import numpy as np import json import os from pyecharts import Map import re #路径目录 baseDir = ''#当前目录 staticDir = os.path.join(baseDir,'Static')#静态文件目录 resultDir = os.path.join(baseDir,'Result')#结果文件目录

整理省-pm2.5表格作为观测值矩阵,将省内所有市的pm2.5的算术平均值作为该省的pm2.5，查看前5行

#读取省市字典 with open(staticDir+'/中国省市字典.json','r') as f: china = json.loads(f.read()) #省-pm2.5字典 provincePm = {'province':[],'pm2.5':[]} #读取各市pm2.5，求平均值作为整个省的pm2.5值 for province in china: pm = []#该省所有市的pm2.5列表 for city in china[province]: #文件名没有市这个字 city = city.replace('市','') path = staticDir+'/中国各市2个月的pm2.5/%s.csv'%city if not os.path.exists(path): continue cityPm = pd.read_csv(path,engine="python")['pm25'] pm.extend(cityPm.values.tolist()) #列表不为空才添加 if pm: provincePm['province'].append(province) provincePm['pm2.5'].append(np.mean(pm)) provincePm = pd.DataFrame(provincePm).set_index(['province']) #观测值矩阵 provincePm = provincePm.drop(['西藏藏族自治区','海南省'],axis=0) provincePm.head() pm2.5province北京市37.807018天津市45.245614上海市24.228070重庆市26.368421河南省37.522478

构建空间权重矩阵，该矩阵反应了地区间空间联系，一般距离越近，权重越大，但是对角线为0 这里要调用calDistance函数计算两地距离，该函数的详情可以查看这篇博客程序在得到距离后，取其负3次幂作为两地的权重，查看矩阵前5行5列

#构建空间权重矩阵 if not os.path.exists(resultDir+'/距离矩阵.xlsx'): spaceMatrix = pd.DataFrame({},index=provincePm.index,columns=provincePm.index) for province1 in spaceMatrix.index: for province2 in spaceMatrix.columns: if not np.isnan(spaceMatrix.loc[province1,province2]): continue #两地距离 try: distance = calDistance(province1,province2) spaceMatrix.loc[province1,province2] = distance spaceMatrix.loc[province2,province1] = distance except: spaceMatrix.loc[province1,province2] = np.nan spaceMatrix.loc[province2,province1] = np.nan continue #地点相同，距离取无穷大，不然后面的倒数会报错 if not distance: spaceMatrix.loc[province1,province2] = 1e+10 #删除空行和列 spaceMatrix = spaceMatrix.drop(['西藏藏族自治区','海南省'],axis=0) spaceMatrix = spaceMatrix.drop(['西藏藏族自治区','海南省'],axis=1) spaceMatrix.to_excel(resultDir+'/距离矩阵.xlsx') else: spaceMatrix = pd.read_excel(resultDir+'/距离矩阵.xlsx',index_col=0) #取反距离，次幂可去1~3 spaceMatrix = spaceMatrix**-3 spaceMatrix.iloc[0:5,0:5] 北京市天津市上海市重庆市河南省province北京市1.000000e-306.807960e-168.255561e-193.228904e-194.260133e-18天津市6.807960e-161.000000e-301.150491e-183.355669e-195.412329e-18上海市8.255561e-191.150491e-181.000000e-303.316709e-191.808123e-18重庆市3.228904e-193.355669e-193.316709e-191.000000e-301.415375e-18河南省4.260133e-185.412329e-181.808123e-181.415375e-181.000000e-30

现在有了省-PM2.5表格作为观测值矩阵，空间权重矩阵，输入到moranI函数中可以得到输出结论并自动绘制莫兰散点图，其中moranI函数位置和代码在下文给出。莫兰散点图，横纵坐标一般是进行分析的数据的离差值，莫兰散点图经常用来研究局部空间不稳定性，其四个象限分别对应于区域单元与其邻居之间四种类型的局部空间联系形式。

#计算莫兰指数 result = moranI(spaceMatrix,provincePm) print('全局莫兰指数为',np.round(result['I']['value'],2),'说明地区之间的pm2.5正相关') print('全局莫兰检验数为',np.round(result['ZI_N']['value'],2),'>1.96，在0.05的水平上显著，认为pm2.5存在空间相关性') 全局莫兰指数为 0.49 说明地区之间的pm2.5正相关全局莫兰检验数为 3.75 >1.96，在0.05的水平上显著，认为pm2.5存在空间相关性

上文中用到的moranI函数位置：Module/MoranI.py，该代码可以直接运行进行测试

# -*- coding: utf-8 -*- import matplotlib.pyplot as plt import numpy as np import os def moranI(W,X): ''' W:空间权重矩阵 X:观测值矩阵归一化空间权重矩阵后进行moran检验，实例https://bbs.pinggu.org/thread-3568074-1-1.html ''' W = np.array(W) X = np.array(X) X = X.reshape(1,-1) W = W/W.sum(axis=1)#归一化 n = W.shape[0]#空间单元数 Z = X - X.mean()#离差阵 S0 = W.sum() S1 = 0 for i in range(n): for j in range(n): S1 += 0.5*(W[i,j]+W[j,i])**2 S2 = 0 for i in range(n): S2 += (W[i,:].sum()+W[:,i].sum())**2 #全局moran指数 I = np.dot(Z,W) I = np.dot(I,Z.T) I = n/S0*I/np.dot(Z,Z.T) #在正太分布假设下的检验数 EI_N = -1/(n-1) VARI_N = (n**2*S1-n*S2+3*S0**2)/(S0**2*(n**2-1))-EI_N**2 ZI_N = (I-EI_N)/(VARI_N**0.5) #在随机分布假设下检验数 EI_R = -1/(n-1) b2 = 0 for i in range(n): b2 += n*Z[0,i]**4 b2 = b2/((Z*Z).sum()**2) VARI_R = n*((n**2-3*n+3)*S1-n*S2+3*S0**2)-b2*((n**2-n)*S1-2*n*S2+6*S0**2) VARI_R = VARI_R/(S0**2*(n-1)*(n-2)*(n-3))-EI_R**2 ZI_R = (I-EI_R)/(VARI_R**0.5) #计算局部moran指数 Ii = list() for i in range(n): Ii_ = n*Z[0,i] Ii__ = 0 for j in range(n): Ii__ += W[i,j]*Z[0,j] Ii_ = Ii_*Ii__/((Z*Z).sum()) Ii.append(Ii_) Ii = np.array(Ii) #局部检验数 ZIi = list() EIi = Ii.mean() VARIi = Ii.var() for i in range(n): ZIi_ = (Ii[i]-EIi)/(VARIi**0.5) ZIi.append(ZIi_) ZIi = np.array(ZIi) #moran散点图 #用来正常显示中文标签 plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示负号 plt.rcParams['axes.unicode_minus']=False fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.spines['top'].set_color('none') ax.spines['right'].set_color('none') ax.xaxis.set_ticks_position('bottom') ax.spines['bottom'].set_position(('data', 0)) ax.yaxis.set_ticks_position('left') ax.spines['left'].set_position(('data', 0)) WZ = np.dot(Z,W) ax.scatter(Z,WZ,c='k') x1 = range(int(Z.min()),int(Z.max()+1)) y1 = range(int(Z.min()),int(Z.max()+1)) ax.plot(x1,y1,'k--',label='x=y') x2 = list(range(int(Z.min()),int(Z.max()+1))) y2 = np.array(x2)*I[0][0] ax.plot(x2,y2,'k-',label='I*x=y') ax.legend(loc='upper right') imgPath = os.path.join(os.getcwd(),'莫兰散点图.png') ax.set_title('莫兰散点图') fig.savefig(imgPath) return { 'I':{'value':I[0,0],'desc':'全局moran指数'}, 'ZI_N':{'value':ZI_N[0,0],'desc':'正太分布假设下检验数'}, 'ZI_R':{'value':ZI_R[0,0],'desc':'随机分布假设下检验数'}, 'Ii':{'value':Ii,'desc':'局部moran指数'}, 'ZIi':{'value':ZIi,'desc':'局部检验数'}, 'img':{'path':imgPath,'desc':'莫兰散点图路径'} } if __name__ == "__main__": w = [ [0,1,1,0,0], [1,0,1,1,0], [1,1,0,1,0], [0,1,1,0,1], [0,0,0,1,0] ] w = np.array(w) x = [ [8,6,6,3,2] ] x = np.array(x) print(moranI(w,x))

最后根据上文得到的局部莫兰指数，绘制中国各省的热力图

#热力图 provice = list(re.sub(r'市|省|.族自治区|维吾尔自治区|自治区','',index) for index in spaceMatrix.index) values = list(result['Ii']['value']) map = Map("pm2.5莫兰指数热力图", '', width=1200, height=600) map.add("", provice, values, visual_range=[-1, 1], maptype='china', is_visualmap=True, visual_text_color='#000') map.render(path=resultDir+"/pm2.5莫兰指数热力图.html")

【本文地址】

公司简介

联系我们